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Abstract 

The distribution of solutions of the Thouless- Anderson-Palmer equation is studied by exten- 
sive numerical experiments for fully connected 3-body interaction Ising spin glass models in a 
level of annealed calculation. A recent study predicted that when the equilibrium state of the 
system is characterized by one-step replica symmetry breaking, the distribution is described by 
a Becchi-Rouet-Stora-Tyutin (BRST) supersymmetric solution in the relatively low free energy 
region, whereas the BRST supersymmetry is broken for higher values of free energy (Crisanti et 
al., Phys. Rev. B 71 (2005) 094202). Our experiments qualitatively reproduce the discriminative 
behavior of macroscopic variables predicted by the theoretical assessment. 

PACS numbers: 05.50+q, 64.60. Cn, 75.10.Nr 
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I. INTRODUCTION 



The Thouless-Anderson-Palmer (TAP) approach is a representative method for investi- 



gating disordered systems of mean field type 



l( . Unlike the replica method, which is another 



method employed in the research of disordered systems for examining average behavior, this 
scheme makes it possible to analyze a system realized by a specific sample of quenched dis- 
order. Using an appropriate approximation of mean field type, the TAP approach yields a 
set of coupled nonlinear equations, which is termed the TAP equation. Solutions of the TAP 
equation provide approximate averages of state variables for the specific quenched disorder, 
which are expected to provide accurate values in the limit of the infinite system size when 
adequate treatment is provided. 

For many systems, the TAP equation can have multiple solutions, the number of which 
is expected to grow exponentially with respect to the systems size when frustration becomes 
sufficiently strong in low temperatures. This indicates that exponentially many metastable 
states can exist for an approximate free energy, referred to as the TAP free energy, and is 
considered to be related to the origin of the replica symmetry breaking (RSB). Therefore, 
analytically evaluating the number of metastable states is important. 

More than a quarter of a century has passed since the first analytical estimate of the num- 
ber of metastable states was performed by Bray and Moore (BM) for the Sherrington- 
Kirkpatrick (SK) model [3(. However, with the recent introduction of the concept of Becchi- 
Rouet-Stora-Tyutin (BRST) supersymmetry (SUSY) from quantum field theory, such as- 
sessments have attracted a great deal of attention [4||. This concept shed new light on the 
solution of the BRST SUSY type, which previously attracted little attention, and on the 
solution of BM, for which the BRST SUSY is broken. A recent theoretical study indicated 
that the BRST SUSY solution is never relevant for the SK model in the level of annealed 
assessment j^J, as verified by extensive numerical experiments Unfortunately, solving 



cper 



the TAP equation is technically difficult in general |2(. Therefore, as far as the authors 
know, numerical validation of the theoretical results has yet been not sufficiently performed, 
particularly for p-body interaction systems, which are considered to exhibit an RSB that 
differs from that of the SK model. 

Thus, we herein investigate the p-body Ising spin models of mean field type numerically, 
in particular, for p = 3. More precisely, we perform extensive numerical experiments for 
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TAP equations of 3-body systems and compare the results with the theoretical predictions 
for the level of the annealed assessment, which has been obtained in a previous study 
The previous investigation utilizing the BRST SUSY concept indicates that the distribution 
of the metastable states of p(> 3) systems in the RSB regime is described by the BRST 
SUSY and BRST SUSY breaking solutions for low and high free energy regions, respectively, 
whereas only the BRST SUSY breaking solution is relevant in the case of p = 2. The main 
objective of the present study is to numerically verify the relevance of the BRST SUSY 
solution of an RSB phase for p > 3 systems. 

The present paper is organized as follows. In the next section, we briefly review how the 
BRST SUSY concept is introduced in counting the number of TAP metastable states and 
introduce the model considered herein. In Section 3, the scheme of the numerical experiments 
is explained and the obtained results are compared with the analytical predictions. The final 
section is devoted to a summary. 



II. TAP COMPLEXITY AND ANALYTICAL EVALUATION 



A. Model and TAP equation 

Although the numerical study is performed only for 3-body systems, we describe the 
general model for p-body systems because the results obtained from the experiments are 
considered to be qualitatively relevant for general cases in which p > 3. The p-body spin 
glass model is defined by the following Hamiltonian: 

H{S) = — Jiiiz-ipSiiSii ■ ■ ■ Si p , (1) 

(iii2—i P ) 

where S = (Si) = {+1,-1}^ and J^...* is the quenched random coupling drawn from 
a Gaussian distribution of zero mean and the variance p\/(2N p ~ 1 ) for each ordered set of 
indices (ziZ2 . . .i p ). The conventional procedure using the TAP approach provides the TAP 
equation, as follows: 

tanh-^m;) = ^ ^ -A.,,,,...,, • • • ™ jp - m^-^ — ^ (1 - q)q p ~ 2 , (2) 

which represents the extremum condition of the TAP free energy 

1 1^/1 + rrii 1 + to; 1 - mi l-mA 

— F T Ap(m) = — > In In 

N v J N ^ \ 2 2 2 2/ 



-jj J iH2...i p m il m i2 ...m ip - —((p-l)q p -pq p 1 + 1) (3) 

(iii2...i p ) 

where q = A -1 $2 i=1 m ?- 



B. BRST SUSY in Distribution of TAP Metastable States 



We here briefly summarize how the BRST SUSY concept is related to the TAP solutions 



. For a specific sample of {Ji t i 2 ...i p }, the number 



following the argument presented in Ref. 
of solutions that have a value of free energy density / can be evaluated as 

TV 

M = ^n^-O^ApM-A/) 



a i=l 

N 



= I drnf[6(d l F TAP (rn))\det(d l d J F TAP (rn))\5(F TAP (rn) - Nf), (4) 
J i=i 

where d{ = d/drrii and a denotes the index of metastable states. In order to proceed further, 
we replace | det(<9j<9jF T Ap(m))| with det (did jF TA p(m)), expecting that the Hessian of most 
solutions contributing to the distribution is positive definite, and represent the determinant 
using fermionic variables {?pi,?pi} as 

det(d^F TAP (m)) = J Vipltyexp (^^jd^F^rn^j . (5) 

Using Eq. (J2J) in the last delta function on F T Ap(m) of Eq. (jlj) to eliminate the J-dependence 
in F^ A p(m), which considerably reduces the work involved in the calculation but is valid 
only on the TAP solution, another expression of p(f) is obtained, as follows: 

p(f) = J VxVrnV^V*pe^ m ' x ^^\ (6) 

S(m,x,il>,il)) = ^(xi mi)diF TAP (m) 

i ^ 

+ ^^^F TAP (m) + u (F TAP (m) - Nf) , (7) 

ij 

where T>a = prefactor x FJ^ dai and u is determined so as to extremize the integral. 

The annealed distribution, which is considered in the present paper, is obtained by aver- 
aging Eq. © directly with respect to {■7i 1 j 2 ...i_}- The obtained result can be characterized 
by the annealed TAP complexity 

£(/) = ^hxp(f) 



= BS.. {~ fuP ~ (1 " q)(B +A) - qX + + k/ 

where p(/) is the average of Eq. (jlj) with respect to { Jj 1 i 2 ...i J ,}, Extr a {---} denotes the 
extremization of • • • with respect to a, \i — /3 2 p/2, 



dm ( 1 



m 2 



+ £? exp 



(tanh (m) — Am} 2 



2/iq 



p-i 



+ Am + u(3f(m) 



(9) 



and 



1 p — 1 

(3f(m) = — ln(l — m 2 ) — In 2 H mtanh _1 (m) 

2 p 

_^[i + ( p -2)(r l -(p-i)g p ]- (io) 

Note that Eq. (jJJ) is invariant under the generalized BRST transformation 5m« = e0*; 
5xi = —eu(p — l)/pipi] 8ip i = —exf, Sipi = 0, where e is an infinitesimal fermionic con- 
stant. Such a property is referred to as the BRST SUSY. This guarantees that for any 
functions of {m, x, if), ip}, 0(m, x, ip, if}), the average over the weight defined by action 0, 
(0(m,x, ■«/>,?/>)), is also invariant under the BRST transformation, which can be used to 
yield various identities. In particular, two identities 

= - ( m i%i) + - ( m l) > (11) 
u {ipi^Pi) = (x 2 ) - y (rriiXi) + ^ (m 2 ) , (12) 

which are obtained by setting O = m i ip i and O = x^ i , respectively, which is useful for 
solving the rather complicated extremization problem of Eq. (|5|). Specifically, by imposing 
two conditions A = —f3 2 {p — l)uq p ~ 1 /2 and A = (3 2 {p — l) 2 M 2 g p_1 /(4p), which represent two 
relations that are obtained by averaging the BRST SUSY identities of Eqs. (fTTj) and (JT2*Jl 
over {Ji 1 i 2 ...i p } with appropriate weights, on Eq. (JBJ and setting B = following BM, a 
considerably simplified expression of the TAP complexity is obtained |4( , as follows: 



Sbrst(/) = Extr {/ftx[F 1RSB (# ?, -«) - /]} • (13) 



Here, 



(3F lRSB (fc q u x) = ~ In J Dz (2 cosh (Pyj^^ ) - ^j^V A? - ^f-, (14) 
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ution (JT3J) is 
. However, 



in 
6] 



where Dz = dzexp[— z 2 /2]/y / 27r, indicates the one-step RSB (1RSB) variational free energy 
of the current p-body spin glass systems, the extremization of which provides the values of 
self-overlap q\ and Parisi's RSB parameter x in the equilibrium state. 

In addition to Eq. (j!3|) . the extremization problem of Eq. (jHj) yields another solution, 
which does not satisfy the BRST SUSY relations. This solution, which we refer to as Srm(/), 
agrees with that which was previously discovered by the conventional recipe of BM |2j. For 
p = 2, S B m(/) always dominates S B rst(/) and, therefore, the BRST SUSY so 
not physically relevant, which is also supported by the numerical experiments 
for p > 2, in particular, p = 3, Sbrst(/) overcomes £bm(/) in a lower free energy region, 
implying the physical relevance of the present approach. We herein focus primarily on the 
justification of this theoretical prediction by extensive numerical experiments. 

III. NUMERICAL EXPERIMENT 
A. Numerical Scheme 

Computational cost is a major difficulty in numerically solving the TAP equation for 
p > 3. The cost for assessing the right-hand side of Eq. (J2J) increases 0(iV p_1 ) per element, 
which practically reduces the upper limit size of systems that can be numerically examined 
as p grows. Since extensive experiments have only been performed up to N = 80, even 
for the case of p — 2, reducing the computational cost to the greatest extent possible is 
indispensable. 

Solving the TAP equation by natural iteration of Eq. (j2J) is not possible Instead, 
gradient descent schemes of a master function L(m) = YliLi (djFTAvim)) 2 , by which con- 
vergence to a certain state is guaranteed, have frequently been employed [l(J [ill Q ■ 
However, simple gradient methods generally require a number of iterations as the locally 
determined direction of steepest descent does not necessarily approximate well the direction 
to the nearest local minimum, which requires several updates for convergence to the local 



minimum. A solution to this difficulty is t 



used in the investigation of p = 2 systems [ft]. Unfortunately, the use of this scheme is also 
difficult for our purpose because assessment of the Jacobian required in this algorithm has 
a higher computational cost than that at the practically feasible level, even if the efficient 



re use of a quasi Newton method, which was 
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formula proposed by Sherman- Morrison is used jl^j ]. 

In order to avoid such difficulties, we employed a simple discrete algorithm that guarantees 
a reduction in the TAP free energy at each update as follows. The diagonal elements of the 
Hessian of the TAP free energy are provided as follows: 



5? Wm) = -L-s + P2P(P - q)q^ > 0. (15) 

J. / f t Zi 



This indicates that solving Eq. (J2J) with respect to a single variable while maintaining 
the other variables {rrtj-ti} constant yields a unique solution, which can be searched by an 
efficient scheme such as the bisection method. In addition, updating to the solution 
guarantees reduction of the TAP free energy (J3J). This implies that for an arbitrary initial 
state, iteration of the above procedure leads to convergence of a certain local minimum of 
the TAP free energy (jHJ). Hopefully, applying this algorithm to a sufficient number of initial 
states while randomly permuting the order of updates among elements will yield a complete 
set of local minima for a given fixed set of {</i 1 i a ...i p } ■ 

A major drawback of this strategy may be that extrema of other types, namely, saddle 
points and maxima cannot be obtained, although such solutions are taken into account 
in Eq. (@J). If the contribution of saddles and maxima to the complexity are relevant, 
this implies that the proposed strategy is not appropriate for verifying Eq. (@J) or Eq. 
(JHJ). However, recent research into the case of p = 2 revealed that most solutions of the 
TAP equation consist of a pair of a local minimum and a saddle point, for which only one 
direction has a negative eigenvalue This implies that the complexity composed of 

only local minima yields the same result as that evaluated from all extrema, justifying a 
crucial approximation | det(didjF^Ap(m))\ ~ det(didjF^p(m)) a posteriori. We, therefore, 
employ the minimization algorithm with the expectation that this property also holds for 
cases in which p > 3, which should be investigated in the future. 



B. Comparison with Theoretical Predictions 

We performed numerical experiments of p = 3 systems for iV = 20, 30, 40 and 50 at 
T = T s = 0.5, which corresponds to the region of 1RSB in equilibrium analysis. The number 
of samples of {J h i 2 ... ip } was 25800 for N = 20, 17200 for N = 30, 1720 for iV = 40 and 430 for 
N = 50. Unfortunately, it was difficult to efficiently find solutions other than the trivial para- 
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FIG. 1: Annealed complexity vs. free energy density for various values of N. The solid and 
dotted lines represent £brst(/) and Sbm(/), respectively. £brst(/) agrees with Sbm(/) and is 
terminated at critical point f c = —0.800 .... 

magnetic solution m = from randomly chosen initial states, which requires a practically 
infeasible number of trials. In order to raise the efficiency of finding non-trivial solutions, we 



first solved the naive mean field equation tanh 1 i 



mi 



a 



(p-i) 



of T = T s from many initial states, which can efficiently produce a number of non-trivial so- 
lutions. The obtained solutions are employed as initial states for solving Eq. (j2J. For finding 
solutions uniformly, we randomly permuted the order of updates among elements at each 



iteration indexed by t; the convergence was judged by a condition J2iLi 
In order to ensure that our search of solutions is as exhaustive as possible, we repeated the 
search 500 times for N = 20, 500 times for N = 30, 2000 times for N = AO and 5000 times 
for N = 50, respectively, for a given sample of connections {Ji 1 i< 2 .„i p \- All solutions found in 
samples were found at least five times. The Hessian was positive definite for all the solutions, 
which indicates that our simple scheme successfully worked as a minimization algorithm. 

Based on a previous study [6], we numerically assessed the complexity by a formula 
£(/) = In \p(f)M] /N, where p(f) is the probability of finding a solution of the free energy 



m\ — m. 
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FIG. 2: Parisi's 1RSB parameter x vs. N. The crosshairs denote numerical estimates assessed 
from slope at the left terminal of the complexity // for which S(/^;iV) = by quadratic fitting. 
The line denotes the theoretical value for T = T s = 0.5, x = 0.710 .... 



density / among the total set of numerically found solutions and M is the average number 
of solutions for one sample set of connectivities. Theoretical assessment predicts that the 
BRST SUSY solution £brst(/) slightly dominates the BM solution £bm(/) m the region 
of lower free energies / < f c = —0.800 . . ., while Sbm(/) is dominant for / > f c . 

Figures Q show a comparison of complexity between the theoretical predictions and nu- 
merical data. These exhibit reasonable consistency in the lower free energy region, where 
Sbrst(/) is dominant, except for N = 20, while the discrepancy is considerable for larger 
/. Solutions of high free energy have been reported to be difficult to obtain (f| 12, 

HQ, 

which is thought to be the reason why in the present study the complexity of larger / is 
underestimated as N increases. Based on the formal equivalence between the BRST SUSY 
solution of the TAP framework and the 1RSB solution in the replica analysis, the tangent 
slope at the left terminal of the complexity /; = —0.823... for which £brst(/j) = is 
expressed as f3x, where x is the Parisi's 1RSB parameter, the value of which is determined 
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FIG. 3: Plefka's stability parameter vs. free energy density for various values of N. The solid and 
dotted curves represent the values assessed by the BRST SUSY and BM solutions, respectively. 

by the extremization condition of the 1RSB free energy ()14|). Figure El shows that numer- 
ical estimates of x assessed by the quadratic fitting of the complexity is in relatively good 
agreement with the theoretical value of x = 0.710 . . . for T = T s = 0.5 except for iV = 20. 

Tempering local minima while maintaining the relation (J2J) from T = to T s might be 
another scheme for improving the performance of the solution search [ijj]. However, such 
a strategy was found to be not so effective, particularly for finding solutions of lower free 
energies. At T = 0, the complexity of this system is expressed by the BM solution of BRST 
SUSY breaking in the entire region of /. Therefore, it can be expected that solutions of 
T = T s , which are adiabatically linked to those of T = 0, are mostly characterized as BRST 
SUSY breaking unless the properties of the solutions are changed. This may be the reason 
why solutions of lower free energy at T = T s , which are characterized by the BRST SUSY 
complexity, are difficult to find by the tempering scheme. 

It is difficult to observe the transition at / = f c between the BRST SUSY and BM 
solutions by numerically assessing the complexity because the difference between Ebrst(/) 
and Ebm(/) is considerably small, even in the region where £brst(/) dominates £bm(/)- 
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FIG. 4: Self-overlap vs. free energy density for various values of N. The solid and dotted curves 
represent values assessed by the BRST SUSY and BM solutions, respectively. 

However, this is not the case for other quantities. 

Figure El shows a plot of Plefka's stability parameter [r| 

(3 2 p{p-l) 



1 - 



((1-m 2 ) 2 ) 



(16) 



versus free energy density /. Theoretical assessment indicates that x p vanishes to zero at 
f = f c signaling the BRST SUSY-SUSY breaking transition, while it is maintained positive 
at other locations (curves). Although the deviation from the theoretically obtained curves is 
relatively large due to the finite size effect, the numerical data are minimized in the vicinity 
of / = f c , which is qualitatively in accordance with the above mentioned behavior. Similar 
behavior is also observed for self-overlap q (Figure 0J). 



IV. SUMMARY 



In summary, we have numerically examined the properties of the TAP solutions for p- 
body Ising spin glass systems, in particular, for p = 3 cases. A theoretical study reported 
that the complexity of this system is expressed by the BRST SUSY solution in the lower free 
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energy region when the equilibrium state is expressed by the 1RSB solution, while the BM 
solution of BRST SUSY breaking is relevant for the higher free energy region |8j . As another 
study concluded that the BRST SUSY solution is always irrelevant for p = 2 cases |6J , the 
emergence of the BRST SUSY solution in the lower free energy region is a discriminative 
feature of p(> 3)-body systems. 

In order to justify this feature, we have performed extensive numerical experiments em- 
ploying a simple discrete algorithm, which locally minimizes the TAP free energy (J3J). A 
heavy computational cost practically limits the experiments for systems of relatively small 
sizes of less than N = 50, which makes it difficult to quantitatively distinguish the BRST 
SUSY solution from that of BM in complexity curves by numerical methods. However, other 
quantities, such as the Plefka's stability parameter x p and self-overlap q, are expected to ex- 
hibit peculiar behavior as a function of free energy density / as a consequence of the BRST 
SUSY-SUSY breaking transition, which has been qualitatively reproduced by our numerical 
experiments. This may indicate that the TAP solutions of lower free energy of p(> 3)-body 
systems in 1RSB phase, which govern equilibrium properties, are characterized as having 
BRST SUSY, unlike p = 2 systems. 

Application of the current approach to sparsely connected systems is left as a topic for 
future study jll|. 
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